Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 5 February 2008 (MN style file v2.2) 



The impact of mass loss on star cluster formation. I. 
Analytic results 
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ABSTRACT 

We study analytically the disruptive effect of instantaneous gas removal from a cluster 
containing O stars. We setup an iterative calculation based on the stellar velocity 
distribution function to compute the fraction of stars that remain bound once the 
cluster has ejected the gas and is out of equilibrium. We show that the stellar bound 
fraction is a function of the initial cluster distribution function as well as the star 
formation efficiency, e, taken constant throughout the cluster. The case of the Plummer 
sphere is dealt with in greater details. We find for this case that up to ~ 50% of the 
stars may remain bound when e assumes values < i , contrary to expectations derived 
from the virial theorem. The fraction of bound stars is expressed algebraically for 
polytropic distribution functions. 

Key words: Stars: formation; star clusters: formation, evolution 



1 INTRODUCTION 

Most if not all stars are formed in clusters or associations, and therefore this is the predominant mode of star formation 
contributing to the Galactic-field population (see reviews in Grebel & Brandner 2002). To understand how such aggregates 
form and evolve remains a severe challenge however, since no theory exists yet that accounts in detail for the (presumably) 
simpler case of the formation of individual stars. On the observational front, surveys of star-forming regions suggest that the 
mass-fraction of gas used up at birth does not exceed 30% of the total when averaged over typical cluster scales (Lada 1999). 
This low star-formation efficiency (SFE) implies that young stars or clusters of stars should be embedded in gas (see Lada & 
Lada 1991). Yet populous clusters with ages J; a few Myrs (e.g. the Orion Nebula Cluster, R136 in 30 Doradus) are already 
void of gas. The mass-fraction of gas left after the epoch of star formation depends crucially on the star formation history 
of the cluster. The formation of proto-stars by hydrodynamical collapse points to the rapid formation of opaque stellar cores 
followed by time-dependent accretion, with the more massive cores accreting at higher rates (e.g. Boily & Lynden-Bell 1995). 
Stellar masses are then accrued over a period of typically 10 s - 10 6 years (Low & Lynden-Bell 1976; Shu 1977; Wuchterl & 
Klessen 2001). Observational evidence in support of rapid accretion as opposed to a slow star formation process sustained 
by ambipolar diffusion of magnetic field is discussed by Hartmann et al. (2001), who point out that the age spreads between 
cluster members are too small compared with that expected from ambipolar diffusion (see Jones et al. 2001 and references 
therein). Because opaque stellar cores will not all have the same mass (or indeed chemical composition), it is unlikely that 
stars will end their accretion phase at the same universal time. Consequently, the star cluster will not be entirely depleted of 
gas when the more massive stars reach the Hydrogen-burning phase. 



OB stars provide an efficient mechanism by which gas is driven out. Stellar winds from OB stars yield a momentum flux 
T ~ 8 x 10 _3 Mq kmsec _1 yr _1 (ChurchweU 1999). This is sufficient to unbound gas from an 1O 3 M0, 1 pc radius cluster 
in w 10 s years, which is comparable to its crossing time t CI — 2R/a oc l/^/Gp(0). Here p(r) is the mass density and a the 
velocity dispersion. In addition, the ionising radiation heats the gas to 10 4 K, causing an over-pressure and expansion at the 
sound velocity « 10 km/s. Furthermore, the star-formation timescale of ~ 10 6 yrs far exceeds a single cluster crossing time 
t cr : the magnitude of stellar velocities should therefore reflect the depth of the cluster gravitational potential since protostars 
would have time to explore the entire cluster volume. This suggests that a rich cluster is already close to virial equilibrium 
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when OB winds set in. 

The question whether or not the stars will re-establish a bound gas-free equilibrium has bearings on the spatial distri- 
bution and statistics of stellar populations in the Galaxy (Kroupa 2002, Grebel & Brandner 2002). A long-standing goal has 
been to understand in what way gas-ejection and the low SFE inferred from observations set constraints on the star cluster 
properties at birth (see e.g. Lada 1999). Hills (1980) has argued from this standpoint and using the virial theorem that if the 
SFE is less or equal to 50% the entire cluster will dissolve. His classic argument holds for clusters that have achieved virial 
equilibrium and which then undergo a sudden loss of mass. Theoretical and numerical work support this conclusion when gas 
evacuation proceeds rapidly (Tenorio-Tagle et al. 1986; Geyer & Burkert 2001), while the numerical N-body study of Lada et 
al. (1984) already found bound stellar cores regardless of the gas evacuation time-scales. This was also found to be the case 
in recent, larger-N studies (Goodwin 1997; Kroupa, Aarseth & Hurley 2001). Adams (2000) has recently pointed out that if 
the SFE peaks with the central cluster density, a small core of stars will always remain bound, even for an SFE well below 
the reference 50% obtained from the virial. This effectively solves the problem of forming bound clusters but at the expense 
of simplicity, by assuming a radial variation for the SFE, the nature of which is yet poorly constrained. Our goal is to show 
that one may solve for the bound fraction of stars analytically without introducing further complications. 

The numerical N-body results and the analysis of Adams (2000) hint that Hills' (1980) original argument is lacking an 
ingredient to account for bound clusters at low SFE. We note that the derivation based on the virial theorem is a statement 
about the global properties of the cluster with regards to the impact of the SFE. A refinement is therefore possible if we 
reformulate the problem from the point of view of the equilibrium stellar distribution function (=DF), F(x,v), as follows : 
First, select stars by binding energy (E) to identify bound (E < 0) and candidate escaping (E > 0) stars. The DF may then 
be modified to take account of the cluster mass loss due to unbound gas and stars. As a second step, we setup an iterative 
scheme to identify bound stars self-consistently once gas and escapers are subtracted from the gravitational potential. This 
scheme allows us to solve analytically for the surviving fraction of bound stars. We show that the result is function of both SFE 
and the chosen DF by contrasting results obtained for polytropes, and in particular the Plummer sphere. The time-evolution 
of clusters, and applications to other DF's, are deferred to a companion paper (Paper II). 



2 HILLS' PROBLEM 

A cluster of stars forms converting a fraction e of the gas into stars in the process. Stellar winds or a supernova event blow/s 
out the remainder on a time-scale r < t cr . What fraction of the stars remain to form a bound cluster? Since the star formation 
epoch will have lasted as long as or longer than a cluster dynamical time t CI , the system as a whole is close to virial equilibrium 
before gas is expelled. In the following, we consider the limit r — > 0, which defines an initial value problem for an open cluster's 
mass profile and velocity field. 

Hills (1980) argued that equilibrium, self-gravitating stellar systems would dissolve if more than half the mass is lost 
instantaneously: stars then preserve their kinetic energy established under a deeper potential well, hence may escape if their 
binding energy E becomes positive. We write the total energy 

E = -^f + \M{J)<Q, (1) 

where M, R are the mass and gravitational radius of a spherical distribution of mean square velocity {v 2 ) . Before gas-expulsion 
the total mass in gas and stars is M = Mi n it = M gaB + M*, while after the gas is expelled M = M t but the velocity dispersion, 
{v 2 } — GAfinit/ 'R, remains unchanged. The new, total system energy is now 

GM 2 1, , , 2x GM 2 

where the last equality follows from applying the virial theorem to the stellar system of mass after it has reached a new 
equilibrium, of energy E* but with radius R* > R. We may rearrange (^[) to obtain 



R* _ 1 M+ _ 1 M-M gas 



7? 2M*-§M 2±M-M g . 



(3) 



The system radius i?* — * oo when A'/ ga s = M/2 and hence the stellar system has zero or positive total binding energy when 
the mass is reduced by 50% or more (e = M*/M < 1/2). 

Hills' derivation applies to any mass distribution. It does not, however, say anything specific about the evolution of the 
mass profile except to relate the final virial radius R+ to the initial system radius R. Our intuition is that a fraction of the 
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stars with larger velocities will escape before we reach the reference 50% gas mass loss and hence the residual stellar potential 
should not be evaluated from the global stellar mass fraction alone. In other words, the stellar velocity DF should determine 
which stellar orbit remains bound according to its binding energy: this is confirmed when considering orbits in a harmonic 
potential as shown in Appendix A. The example of harmonic motion will serve us later in interpreting results obtained for 
systems with realistic mass profiles and DF's. Below we introduce a scheme to pin down more accurately the relation of bound 
mass to star formation efficiency, e. 



3 DISTRIBUTION FUNCTIONS 
3.1 Iterative scheme 

Consider an isotropic mixture of stars and gas within a spherical volume. The stellar DF F[r, v) is then a function of energy 
alone. As a simplifying assumption we take the SFE to be independent of position in the cluster. Hence the SFE parameter 
e = constant. Thus gas and stars are initially mixed in the same proportion throughout, however note that e does not fix the 
profiling of stellar masses with radius, rather the total mass of gas made into stars : 

V\ (number of stars in mass bin i) x (stellar mass mi) 

— ! = constant 

at any radius r. Both gas and stars are distributed according to the same DF. The total system mass is then 

M — M* + Mgas = e _1 = e _1 ff F(r,v)d s rd 3 v = e~ 1 f F(E)W(E)dE, (4) 

J Jv Je c 

where \E C \ is the maximum binding energy and W the density of states of energy E. The cluster potential (j>(r) follows from 

integrating twice Poisson's equation 

0(r; M, R) = f ^ f* 4nGp(w)w 2 dw = y(r/R) (5) 

Joe u Jo 

where y(x) is a dimensionless function. Keeping R constant while removing the gas instantly so the potential includes stars 
only, we have 

4>(r;M,R)^ </>* = £</>, (6) 
and hence the fraction of stars at radius r which now have positive energy is 



F(r,v) d ru du / f(v)v dv 
\ - - "«.« _ ^"e.* ( 7 \ 

/ F(r,v) d i rv 2 dv / f(v) v 2 dv 
Jo Jo 

with v e the (local) escape velocity computed from cf> (before gas- expulsion ) , and similarly for v e>i , obtained from For 
convenience of notation, we write the stellar velocity distribution function f(v) defined at constant radius r as 



f(v) d s v = F(r,v) d 3 r d 3 v . 

For a given DF and e, we may compute the quantity 1 — A e at all radii. Note that should the SFE be a function of 
radius, the re-normalisation (^j) leading to <j!>* would not apply, however (^) may still be computed if the stars' potential is 
given and v e known from the initial gas + stars mixture. The key step is to adjust the potential itself, since potential and 
velocity field do not match any longer. This is normally computed numerically as an initial value problem, for instance using 
N-body integration which takes into account all star-star interactions and the redistribution of energy between them in the 
time-dependent potential. We take a different approach. If we thought that stars acquire only little extra energy during the 
time that they escape, then the fraction of stars remaining might be computed as follows. 

Since the fraction (jjj) of positive-energy stars is known at all r, we re-compute the gravitational potential counting only 
stars with E < at each radius. Neglecting dynamical evolution, the cluster radius R and structure function y(r/R) are 
unchanged and hence the potential can be re-computed by integration, from oo, inwards. Once the new potential is known, a 
fraction of the remaining stars will again be unbound by virtue of the stars lost during the previous iteration. We therefore 
repeat the procedure until finally the cluster mass converges to a finite quantity, in which case no more stars escape and the 
original distribution function is depleted from all escapers in a self-consistent manner. 

We consider the case of scale-free DF in detail, then turn to more physically-motivated cases in the following sections. 
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3.2 Example : scale free DF's 

The case of a power-law isotropic velocity distribution, 

f(v) d 3 v oc v l3 ~ 2 d 3 v oc v 13 dv , 

where (3 is constant, provides a tractable starting point. The DF is truncated at the local escape velocity and we take the 
potential (j>(r) to be a power of the radius: the system is therefore scale-free. The one-dimensional velocity dispersion obtained 
from Jeans' equation which for power-law density and potential yields 

1 f°° 

° 2 (r) = — ~, — r / p(x)V6(x) dx = -2K Mr) = K v 2 (8) 
P[r) Jr 

at each radius for some constant of proportionality, K. Substituting p(r) — > ep(r) in (||), we find the new potential from (^|) 
and escape velocity 



«e,*(r) = V^-MO = e 1/2 v e (r). (9) 



Equation (Q) becomes (with u = v/v e 



2 f 1 

f(v)v dv / f{ u )udu 

1 " Ae = ^ = ^ = 1 - £ (,3+1)/2 , (10) 

^0 Jo 

where we took /3 > — 1. Thus A e = e( /3+1 ' / ' 2 is independent of radius. Repeating the procedure to take account of the 
positive-energy stars, we substitute e — > A e • e, etc, so that after j iterations 

with A e [0] = 1. Assuming convergence of the series with j — > oo we may solve for A e to obtain 

l + (3 

x. = A-P. (11) 

(See Boily & Kroupa 2002 for an alternative derivation.) In terms of the stellar mass fraction, the predicted mass of bound 
stars, M', is then simply 

Ml = \ e eM = \ e M* (12) 

(and similarly for the stellar density p*[r] wrt p[r\). The eigenvalue X e in (|l^) converges to a non-zero (positive) value as 
j ^ oo only for f3 < 1, since e < 1. All DF's with /3 > 1 lead to cluster disruption, because the high-velocity range of the DF 
is too densely populated, leading to catastrophic stellar loss after the expulsion of any amount of gas. 

Fig. 1 graphs the fraction of bound stars M^/M* obtained from (^) as function of the SFE e for five values of the power 
index f3. All curves meet at e = 1 and M| /M* = 1. Note how a flat DF (f3 = 0) predicts a linear decline of bound star fraction. 

The volume density p(r) is recovered from integrating F(r,v) over all velocities at constant radius r : 

p{r)d 3 r= [ ' F(r,v) d 3 r d 3 v = f " f(v) 4tt v 2 dv oc / \ ? dv = = (~2cj>[r]) il3+1)/2 = (-20)' 1 (13) 

Jo Jo Jo p + 1 



where the last relation emphasises the similarity with polytropes of index n, to which we return. Since ( |12| ) was derived 
assuming a power-law potential, it follows from ([l3|) that p(r) is also a power of the radius. Power-law density profiles may 
be of application to young clusters such as NGC2282 (Horner et al. 1997). Letting p(r) oc r~ a , we deduce from (^|) that 
r) oc r 2 ~ a ; a centrally condensed cluster requires a > 0. The power indices a and (3 are then related to each other using 



(E 

a = 2. 

(3-1 

and hence a > for /? > 1 or (3 < —1. We recall that /3 > — 1 from (0), hence we must choose (3 > 1 for a self-consistent 
solution. NGC2282 shows an axially-averaged surface density E oc R~ x which can be obtained from projecting a spherical 
cluster of volume density p oc r -2 , or a = 2 (singular isothermal sphere). The index /3 = (a + 2)/(q — 2) — *■ oo when a = 2 
which is unrealistic. Therefore no self-consistent solution of this sort exists with /3 > — 1 applicable to NGC2282. The large 
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value of (3 obtained for this case is at best indicative of rapid dissolution for any e < 1. 

The example of power-law DF's suggests a strong dependence of the effective SFE when cluster disruption occurs, and 
the shape of the DF itself. Naive application to NGC2282 of a power-law DF based on the Jeans equations of equilibria leads 
to a nonphysical solution. In the following we apply our iterative scheme to other known DF's. 



4 PLUMMER MODEL AND THE FAMILY OF POLYTROPES 

We wish to extend our basic result (jll]) to a class of DF's with properties similar to those of observed star clusters. The 
power-law relation of the scale-free DF's suggests that we consider the family of polytropes defined as 

F{E) = {-E) n ~i (14) 

where the power index n > 1/2 may take non-integer values (BT+87). The case n — 5 has been studied in detail by Plummer 
(1911). Then 

F(E) = {-Ef 2 = (|#r)|- i. 2 ) 7/2 . (15) 
The total system mass of a Plummer sphere is finite but infinite in extent; the density profile 

p(r) = Sf [i+r 2 /*r 5/ WM = (^w as) 

where M is the total system mass, R p the Plummer radius and both the velocity dispersion and density maximise at the centre 
(see e.g. Spitzer 1987). At large distances, the density profile p(r) tx r~ 5 decreases like a power of the radius. Remembering 
(|), the Plummer velocity DF f(v) oc (1 — [v /v,,] 2 ) 7 ^ 2 . Inserting this into (pj) and with 

. , , 1210 2104 2 1488 3 384 4 

be El eH e e -\ e 

yy ' 105 105 105 105 

we find 

/ (1 — u 2 ) 7 ^ 2 u 2 du 

A e = % = 1 - - (sin" 1 e 1 ' 2 - p(e) e 1 ' 2 VT^l) . (17) 

(1 - u 2 ) 7 ' 2 u 2 du 

o 

We note that, as for the case of power-laws, the solution is independent of radius. Therefore the potential 0* and escape velocity 
v e follow from (|^) . The fraction of stars labeled for escape are removed from the potential and the integral re-evaluated with 
the substitution of e 1 ^ 2 by (eA e ) 1,/2 as the lower bound of integration in (|l7]). This follows from the same transformation 
leading to ([]). We may repeat the procedure until S(l — A e ) — > on successive iterations so no additional stars are lost. The 
converged fraction A e is then solution of 

A e = - (sin-^A.e) 1 / 2 - p(A e e) (A e e) 1/2 Vl - A.e) . (18) 
This may be expressed in parameterised form. Defining 

x = A e e ; iC[0,l], (19) 
we have the solution \ e {x) from (|l^). The full solution follows from solving for e as a function of x, 

X e [x) 

which is known for all x in the allowed range. The relation A e (x) to e(x, A e ) saves us from iterating on (^|), which otherwise 
yields the same result. In Appendix B we extend the solution ( |l 8[ ) to configurations out of dynamical equilibrium. 

4.1 Analytic results for Plummer spheres 

From the definition (^|) we must have e = A e = 1 when x = 1. Inserting (^) into ( po| ) we find 



lim e(x) oc x 

x^0 



-1/2 
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while X e (x — > 0) = 0. Since e(l) = 1 and e(x < 1) < 1 for a physically meaningful solution, it must reach a minimum 
somewhere in the interval x : [0, 1]. We find d e/ dx = e(x) — from basic calculus at 



x = 0.225259 ; minimum of e = 0.442201 

when A e = 0.509404. Solutions exist only for an SFE in the range e : [0.442201, 1]. The result @ is displayed on Fig. 2 
where we graph the same quantity A e as in Fig. 1 but for a Plummer DF. The curves are best understood if one starts at 
(e, A e ) = (1, 1), upper right-hand corner of the plot, shifting to the left as the SFE e is reduced. The most striking feature is the 
sudden drop in bound fraction A e for an SFE « 44.22% (solid curve, Fig. 2). The lower branch of the parameter solution X e (x) 
is shown as dash but corresponds to no physical state of the system. For comparison, we have also displayed as open circle 
the results of a series of N-body calculations, which show remarkable agreement with (|l8|). Details of the N-body calculations 
will be found in Paper II. 

Around the critical value, the predicted bound fraction varies from A e ~ 0.63 upward of an SFE e = 0.443, to for an 
SFE < 0.442. For e — 0.50, we compute a fraction of bound stars A e ~ 86%. This counts as bound nearly all stars initially in 
the system, when we expected complete dissolution from ([|). 

To illustrate the significance of repeated iterations on (0), we also plot on the figure the result when no iterations are 
performed and the solution A e (e) taken directly from (^) (dotted line, Fig. 2). This graphs a selection of stars by energy but 
without taking into account the shallower gravitational potential as a result of stars leaving the system. As the bound fraction 
remains finite for any positive SFE, and no sudden dissolution is predicted, we conclude that feedback on the gravitational 
potential from stars lost changes the character of the solution very significantly. This should not come as a surprise since 
cluster dissolution implies a stream of stars moving outwards. What this shows is that the impact of this stream must be taken 
into account when quantifying the fraction of stars which ultimately will form a bound equilibrium after a phase of dynamical 
evolution. By contrast, Adams (2000) allows e to vary with position but does not iterate to establish final membership. 



4.2 Analytic results for all polytropes 

The results of the previous section may be extended to all polytropes of index n > 1/2 defined by (jl4|). The general solution 
X e (x) may be written for any index n according to the definition (1^). We find 

x s/2 2 F 1 (§,§ -n, &,x) _ r(n + l) 4a; 3 '' 2 /3 3 5 



where T(u) is the T— function, and the hyper-geometric function 2-Fi(fc, I, m, z) assumes polynomial expressions for integer 

values of n (see e.g. Weisstein 1999 for details). We may extend the steps for the Plummer sphere to any polytrope, identifying 

the minimum of e(x) in each case. This is done easily with software such as mathematica or MAPLE. 

Results for x where e = are given in Table [j] for several polytropes of index n ranging from n = § to 40. The minimum 

_ l 

of e is also listed along with the stellar mass fraction A e . The case of n = 3/2 is special. Then e = x 2 everywhere in the 
allowed interval, and so e'(x) < for all x's. Indeed no solutions exist with e(x) < 1 for all indices n in the range i < n < |. 
The n = | polytrope admits e = A e = 1 as a special solution satisfying our constraints on bound mass, however e'(l) 7^ 
does not satisfy the mathematical requirement of a local minimum. 

All polytropes of index n > 5 stretch to infinity, and those with n > 5 all have infinite mass. We note that e decreases 
smoothly with increasing n and in particular e — > as 1/n. The same is true of the mass fraction A e , though for n — 40 
we find A e « 0.40 for a relatively low SFE of e as 5%. The bound fraction A e at low SFE recovered from the polytropic DF 
overlaps with observational estimates if the index n is sufficiently large. However, all polytropes with n < 5 dissolve at a larger 
SFE than the Plummer case (n = 5) . Since n < 5 polytropes show a broad region of near-uniform volume density, these data 
suggest that the core-halo structure plays a role in the survival of clusters, in the sense that more concentrated polytropes 
(larger index n) resist dissolution better. 



5 TIME-EVOLUTION AND OTHER CONSIDERATIONS 

The analytic approach of Section 3 did not aim to relate final equilibrium to initial conditions, as we have completely ignored 
details of the time-evolution. We turn here to basic considerations of cluster dynamics that were left out until now. 

Adams (2000) has recently revisited the problem of forming bound clusters through dynamical evolution. His analytical 
approach leads to the conclusion that in order to form bound clusters, of low SFE averaged over the entire cluster, the SFE 
must peak with the central stellar density. In his problem, therefore, the large concentration of stars results from efficient star 
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formation in a small volume. Our own results would support this view, in as much as we find the more centrally concentrated 
polytropes to dissolves at lower SFE. We did not, however, consider an SFE varying radially with the density. The problem 
would seemingly be solved from the freedom in profiling the SFE with the stellar potential. Tenorio-Tagle et al. (1986) had 
noted that gas density gradients have a strong impact on the outcome, since more bound gas (steep gradients) requires 
more mechanical energy for removal, which demands more time under a constant UV flux: this in effect reduces the ratio of 
dynamical time to gas-evacuation time and shifts the dynamics toward a slow evolution regime, which we briefly outline below. 

In condensed star clusters, the dynamical time t CI at the half-mass radius may be long compared to a finite gas-evacuation 
timescale = r, while being shorter than r at the centre of the cluster. When t CI /r <C 1, the integral over an orbital path I 



is an adiabat so J$ = constant during evolution (Weinberg 1993; BT+87, §3.6). Here M is the cluster mass, E the binding 
energy and v, r the orbital velocity and radius of a star, respectively. The regime tcr/i" <C 1 reduces to an adiabatic, quasi-static 
transformation of the initial mass profile, with mass loss acting as a positive source of energy. Notice that in order for ( |22[ ) 
to apply, no limit is set on the fraction of mass lost through gas expulsion, only that it proceeds slowly: then no stars are 
lost through adiabatic evolution. Indeed under such condition clusters prove much more robust to total disruption (Lada et 
al. 1984; Goodwin 1997). From (^) we deduce t CI oc yj r 3 /M oc r 2 /J^. Since the system expands, it will move out of the 
adiabatic regime t CI /r <C 1 (i.e., t CI w r) if the expansion factor 



We can draw constraints from observations of Scorpius clusters which suggest an expansion factor of up to five between 
clusters of different ages (e.g. Brown 2001). If we take r « 10 years as reference (see the introduction), and inserting a ratio 
= 5 on the left-hand side of (|2^), the central region of a cluster would evolve adiabatically throughout the gas expulsion 
phase if the core t CI is initially on the order of 4 x 10 3 years. For a condensed cluster such that initially t CT > r at the 
half-mass radius, we may assume that t cr oc 1/ yj G(p) > r from that point on. The ratio of central to half-mass densities 
p(0)/pi/ 2 = *cr(ri/ 2 )Acr(0) > T 2 /t 2 CI (0) would then be > (10 5 /4 x 10 3 ) 2 = 625 ~ 10 3 . A Michie-King cluster model (cf. King 
1966) with $/<t 2 > 9 meets this condition, however a Plummer sphere does not. This points to yet greater dependence on the 
DF than was obtained for the family of polytropes. Results obtained for the King family of DF will be presented in Paper II. 

This analysis reinforces our view that instantaneous gas removal provides the most severe criterion for the survival of 
bound cluster, since any realistic situations will be a compromise between rapid gas removal and the adiabatic (t CI <C r) 
regime. Hills (1980) had pointed out that explosive mass loss may yet occur at a time when the cluster has not fully reached 
dynamical equilibrium. If the stars' velocities were small initially, the cluster would contract radially until T/| W| ~ 1 so that 
Q — 2T/\W\ = 2 at most (e.g., Boily, Athanassoula & Kroupa 2002). Two possibilities may occur: in the first one, the gas is 
drawn in with the stars so that the mass ratio of gas to stars is unchanged along the path of the stars. In this case we may 
compute the effect of instantaneous gas-expulsion by removing stars with positive energy and setting Q > 1 in the DF (see 
Appendix B). We find for Plummer spheres that total disruption would occur if the SFE were less than about 48.8% when the 
virial ratio assumes its maximum value of Q = 2. Clearly this has the contrary effect to the one sought. A second possibility 
is that the gas is expelled when the stars have not acquired large velocities and Q < 1, a more likely situation given that a 
cluster initially at rest would spend more time in this dynamical phase as can be deduced from Kepler's third law. We find 
in this case for a Plummer sphere and for Q w 1/4 bound stars at an SFE ~ 24%, as in observations of star forming regions. 
This provides a solution for bound clusters where stellar ages are less than a mean dynamical cluster crossing time and virial 
equilibrium has not yet been established. This view has to be weighted against the likelihood of observing clusters out of 
equilibrium and the reduced time available to form stars. 

6 GENERAL CONCLUSIONS 

We have presented a simple method for determining the fraction of stars that will survive a phase of rapid gas loss following 
the birth of the first heavy stars to form a bound cluster. This approach, based on an iterative selection of stars by binding 
energy, is a modified version of the problem discussed by Hills (1980). The star formation efficiency e is taken as free parameter 
which is held constant throughout the cluster. Some highlights from this study are: 

1) The choice of distribution function determines the fraction of bound stars that remain to form a bound cluster following 
the rapid loss of gas (cf. Eq. [?]) . Taking account of the DF allowed us to estimate more accurately the bound mass than the 
predicted disruption for SFE < 50% based on the virial theorem alone (cf. Eq. ||). 

In Appendix A, we show using a harmonic oscillator model that a large fraction of a self-gravitating stellar system in 




(22) 




(23) 
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virial equilibrium will survive gas expulsion despite a low SFE (< 0.4) if the stellar velocity distribution function favours stars 
with low velocities. It is, therefore, important to measure stellar velocities in a cluster-forming molecular cloud core to obtain 
observational constraints on the shape of the DF. 

2) Analytic solutions were obtained for all polytropes of index n > |. Polytropes with index n ^> 1 yield a bound fraction 
for very low SFE, easily matching the range suggested by observations. 

These findings are backed up by a series of N-body calculations, where the orbital evolution of individual stars is fully 
taken into account: such calculations provide a good match to the solution curve ( |l8| ) for the Plummer sphere (Fig. 2). This 
encourages us to apply (12) to other DF and to consider the dynamics fully in order to identify which, if any, initial conditions 
DF may lead to bound clusters despite very low SFE. We propose to cover these topics in greater details in a companion 
paper (Paper II). 
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APPENDIX A: BOUND ORBITS IN A HARMONIC POTENTIAL 

To see how stars on different orbits may escape at different rates, as argued in Section 2, consider a self-gravitating homogeneous 
mass distribution of density p. The stellar orbital time oc t CI is then independent of position. Stars are often distributed 
uniformly in space over a small volume so this situation has bearing on actual systems. Our aim is to relate the membership 
of a star to its orbital parameters prior to gas expulsion. 

The equation of motion for this case reduces to that of a harmonic oscillator which in one dimension reads 

r(t) = V0(r) = GAf (< r) = -to 2 r(t) (24) 



where we set u = 4nGp/3. Integrating (|24j) we have 

r(f) =r a cos(wi + o ) (25) 

where the amplitude r and angle o are set by the initial conditions. When the homogeneous region is cut off at radius r = R, 
the potential everywhere in space is 

3 GM uj 2 2 

^ = ~2~ + ~ r r ~ R > (26) 

and 4>(r) = —GM/r for r > R. The derivative of <f> is discontinuous at r — R but this has no bearing on the argument 
presented here since we consider motion with r(t) < R only. We imagine an ensemble of stars sharing the same orbit of 
amplitude r a but otherwise phase-mixed. The energy per unit mass for that orbit 

£ = <t>(r) + ±v 2 (r) = <j,(ro) 

is modified at time t — t when a fraction of the mass is lost and only stars remain: </>* (r) = e<j>{r) is now the gravitational 
potential. The new energy per unit mass at radius r 

£, = Mr) + \ v 2 (r) = e0(r) +£- 4>{r) 
is positive if v 2 > — 20, (r) = -~2e<f>(r) which is possible provided (regrouping terms and remembering [Efjj ) 

2 2. 2 2 ^ 3 GM 2„2 fn n\ 

uj r a + e u> r >e — - — = e ,iuj R . (27) 
R 

Since the original orbit (^) has r(t) < r a , we find from ( p7| ) a minimal amplitude r for unbound orbits as function of the 
SFE e : 



'" < ^ < 1, (28) 



1 + e ~ R 

where the upper limit follows from our definition of R. Equation (^8[) states that all stars on orbits with amplitude r 
sufficiently small will remain bound, while at least a fraction of those on large- amplitude orbits will escape. Otherwise said, 
stars on low-binding energy orbits, ie which achieve high velocities at some point on their orbit, are more likely to unbind 
following gas expulsion. All orbits with r <C R would lose none or a smaller fraction of stars than those that visit the edge 
of the system. 



APPENDIX B: DEVIATIONS FROM DYNAMICAL EQUILIBRIUM 

The virial ratio Q = 2T/\W\ = 1 for a system in dynamical equilibrium. In practice the DF may not satisfy this condition 
exactly and this requires modifications to (|f7|). We consider the case of the Plummer sphere only but a similar treatment 
would apply to all polytropes. 
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Substituting absolute values for E in ( |15| ) allows to treat sub- as well as super-virial dynamics. Then the integrals (|1 

read 



/ (\l-u 2 \) 7/2 u 2 du 



, ... , HQ) + G(Q) 

( 1 1 — u I ) u du 



1 - A « = = 1 - -r^ , H77^ (29) 



where the function T is defined by 



with p(x) defined by (|17|), and 



.F(a;) = ^sin 1 x 1 ^ 2 — p(x) x 1 ^ 2 Vl — x\ (30) 



S(0 > i) = ggg [x/QV'Q - Tp(0) + iog(>/Q + ^/Q - 1)] (31) 

and Q(Q < 1) = 0. Clearly the numerator in (^) evaluates to zero if e > Q when the solution A e = 1. Since the SFE e 
is bounded to [0, 1], this is possible only for sub-virial conditions (Q < 1). We recover ( Jl7| ) by setting Q = 1 in (|5i]) since 

^(1) = 77/2. 

Sub-virial conditions are relevant to the problem of forming bound clusters (e.g. Klessen & Burkert 2000; Clarke et al. 
2000), with the caveat that sub-virial conditions at the gas-expulsion time require synchronous star formation throughout 
the cluster to within a fraction of the dynamical time, when the spread in colours suggest at least a formation time of a few 
million years (Hartmann et al. 2001). We wish to show that a point of the parameterised solution ( |is| ) and (^) shifts to lower 
e as Q is reduced from 1. Treating Q as the variable in (^) and differencing we have 

d e(zo) _ Xq tt' tn\ ■ d \ e (x ) _ \ e (x ) ^-i (n s ,a \ 

where x , the parameter defined in (^), identifies a point of the solution curve (A e , e) when Q = 1. The derivative 

F'(Q) = ^ y/Q (1 - Q) 7/2 > V Q < 1 . 

It is easy to see from the above equation and (32) that reducing the virial ratio ( dQ < 0) shifts the abscissa e to smaller 
values while increasing the bound fraction A e , for any point x - Therefore 'cooler' initial conditions lead to bound clusters 
for lower SFE e. The gap from finite A e to zero (total dissolution) is more pronounced but occurs at lower values of e as Q is 
reduced. 

Since J-'(Q) is near zero around Q = 1, the net effect of reducing the virial ratio Q from unity is at first marginal and 
only becomes significant once Q < ^. Similar conclusions apply to super-virial (Q > 1) conditions. Total dissolution occurs for 
larger SFE e in these cases. A cluster collapsing from rest would reach Q — 2 at most during violent relaxation. Dissolution 
occurs for an e ~ 48.4%, only marginally larger than the 44.2% value obtained for equilibrium configurations. The solution 
for Q — 2 and e = 1 (SFE of 100%) yields a bound fraction A e « 0.91, such that 9% of the stars escape. This is remarkably 
similar to what is observed in N-body realisations of collapsing cold spheres. 
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Table 1. Critical SFE e and bound fraction A e in parameter form (|l^) for polytropes of index n. 

Comments 



n 


X 


or € 


\ 

A e 


3 


1.000 


1.000 


1.000 


2 


0.528 


0.857 


0.616 


3 


0.470 


0.698 


0.673 


4 


0.296 


0.543 


0.546 


5 


0.225 


0.442 


0.509 


6 


0.182 


0.373 


0.485 


8 


0.131 


0.283 


0.463 


10 


0.102 


0.228 


0.448 


20 


0.044 


0.115 


0.424 


10 


0.024 


0.058 


0.412 


CO 


0.000 


0.000 





SFE close to (g|) 
Plummer sphere 



Isothermal sphere 




Figure 1. Ratio M^/M* of bound to initial stellar masses as function of star formation efficiency for power-law DF (Jl2|). Five values of 
the power index /3 are displayed. Those with fi > are weighed in favour of high-velocity stars: the fraction of bound stars drop more 
sharply with decreasing SFE. 
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Figure 2. Ratio A e = I M* of bound to initial stars as function of star formation efficiency, e. The results for Plummer models derived 
from the iterative scheme jl8^ are shown as the solid line. By contrast, the dotted curve is obtained from (|l7|), when no iterations are 
performed: in this case a finite bound fraction would be expected for any positive SFE. The results of numerical N-body computations 
are shown as open circles. 



